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This paper determines the excess free energy associated with the formation of a spherical cavity 
in a hard sphere fluid. The solvation free energy can be calculated by integration of the structural 
changes induced by inserting the cavity using a number of different exact thermodynamic pathways. 
We consider three such pathways, including a new density route derived here. Structural information 
about the nonuniform hard sphere fluid in the presence of a general external field is given by the 
recently developed hydrostatic linear response (HLR) integral equation. Use of the HLR results 
in the different pathways gives a generally accurate determination of the solvation free energy for 
cavities over a wide range of sizes, from zero to infinity. Results for a related method, the Gaussian 
Field Model, are also discussed. 



I. INTRODUCTION 

The solvation free energy determines how readily a so- 
lute can be dissolved in a given solvent fluid. This plays 
an important role in many chemically and biologically im- 
portant processes, perhaps most notably in hydrophobic 
interactions in water. A significant part of the solvation 
free energy arises from the required expulsion of solvent 
molecules from the region occupied by the harshly re- 
pulsive molecular core of the solute. These very strong 
"excluded volume" interactions can significantly perturb 
the local density around the solute and cause simple ap- 
proaches based on gradient expansions to fail. 

These effects can be seen most clearly in the simple 
model system treated in this paper. We will calculate 
the excess or solvation free energy associated with the 
insertion of a spherical cavity with radius R v into a hard 
sphere fluid, whose molecules have diameter a. By defini- 
tion, the centers of the solvent molecules are completely 
excluded from the region of the cavity, which thus acts 
like a hard core external field. 

This system has many interesting limits. When the ex- 
clusion field or cavity radius R v equals a, then the cavity 
acts like another solvent particle and the solvation free 
energy is directly related to the chemical potential of the 
solvent. As the cavity radius tends to infinity it effec- 
tively turns into a hard wall and the relevant thermody- 
namic quantity is the surface free energy or surface ten- 
sion associated with a hard wall in a hard sphere fluid. A 
cavity with radius R v = a/2 acts like a hard core "point 
solute" of zero diameter. Even shorter-ranged hard core 
fields or "tiny cavities" with R v < a/2 are also of inter- 
est, since the induced structure and solvation free energy 
of a tiny cavity can be calculated exactly. This limit can 
thus serve as a nontrivial check on approximate methods. 

The most commonly used method today for such prob- 
lems is weighted density functional theory (DFT) p|. 
Here one attempts to describe the free energy directly 
as a functional of some kind of smoothed or weighted 
average of the nonuniform and often rapidly varying sin- 



glet density. This has the advantage that the free en- 
ergy is obtained directly and by construction the associ- 
ated fluid structure (obtained by functionally differenti- 
ating the free energy) is consistent with the approximate 
free energy. However the choice of appropriate weighting 
functions is by no means obvious and a number of differ- 
ent and often highly formal schemes have been proposed. 

We focus instead in this paper on making direct use 
of structural information about the nonuniform solvent 
fluid to obtain the solvation free energy. We believe this 
allows physical intuition to play a more central role and 
we can take advantage of the recent development of a 
generally very accurate theory relating the structure of a 
nonuniform hard sphere fluid to the associated external 
field 0. 

As we will see below, the free energy can then be cal- 
culated by integration, starting from an initially known 
state (e.g., the uniform fluid) and determining the free en- 
ergy changes as the solute-solvent interaction (the hard 
core external field) is "turned on" , or alternatively, as 
the density is changed from the initial to the final state. 
There exist many possible routes from the initial to the 
final state, and we will generally refer to them as ther- 
modynamic pathways. If exact results are used for the 
intermediate values of the structure and associated fields, 
then all these different pathways will give the same (ex- 
act) result for the free energy. 

In practice, of course, approximations will have to be 
made and the different pathways will generally yield dif- 
ferent results. This is sometimes referred to as the "ther- 
modynamic inconsistency" of structurally based methods 
. But this can be viewed more positively as giving one 
the freedom to choose particular pathways that could be 
relatively insensitive to the errors that exist in the struc- 
tural theory, and we will try to use this flexibility to 
obtain the most accurate results. Moreover, there is an 
inherent smoothing of the structural information in the 
integration used to obtain the free energy. The differ- 
ences in free energy predicted by different pathways will 
also give us some indication about the overall quality of 
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the theory. 

This approach generally requires the density profiles 
and associated fields of all the intermediate states along 
the various pathways, and thus a fast and accurate 
method for determining these quantities is crucial for 
computation efficiency. We will use here the generally 
accurate hydrostatic linear response (HLR) equation 2] 
proposed by Katsov and Weeks. A different physically 
motivated derivation of the HLR equation is given be- 
low. 

We will also examine the alternative free energy pre- 
dictions that arise from a theory closely related to the 
HLR equation, the Gaussian field model (GFM) devel- 
oped by Chandler Q . For a solute with a hard core the 
GFM proposes an approximate partition function from 
which the associated density response can be derived. In 
the particular case where a rigid cavity is inserted into a 
hard sphere fluid, the HLR and the GFM approaches turn 
out to make identical predictions for the induced struc- 
ture. Thus structurally based routes to the free energy 
involving only hard core fields will give the same results. 
In addition, one can use the approximate GFM partition 
function to evaluate the solvation free energy directly. 
However, as we will show later, the latter approach tends 
to produce less accurate results. This deficiency shows 
up even more strongly in the tiny cavity limit where the 
structural predictions of the HLR and the GFM are ex- 
act, and several pathways giving the exact free energy 
can be found. This illustrates the advantage of consider- 
ing a variety of thermodynamic pathways that can make 
best use of the available structural information. 



II. DENSITY RESPONSE TO AN EXTERNAL 
FIELD 

A. The HLR equation 

We describe the system using a grand canonical en- 
semble, and thus want to determine the excess grand 
free energy arising from insertion of a spherical cavity 
or hard core external field of radius R v . To derive the 
HLR equation we start with the basic linear response 
equation [4| for a nonuniform hard sphere system in a 
general external field 0(r), with chemical potential p B , 
inverse temperature (3 = (fc B T) _1 and associated density 
p(r; p B , [(f)]) = p(r): 

-W(n) = ydrax-^ri.rajWJ^ra). (1) 

This relates small perturbations in the density and field 
through the (inverse) linear response function 

X -1 (ri,r 2 ; [p]) = <5(n -r 2 )/p(r x ) -c(n, r 2 ; [p]). (2) 

Here c(ri , r 2 ; [p] ) is the direct correlation function of the 
nonuniform hard sphere system. The notation [p] indi- 
cates that these correlation functions are nonlocal func- 
tionals of the density p(r). 



Since we want to focus on the effects of the perturbing 
field, we have used the inverse form of linear response 
theory |j| in Eq. QJ, where the field appears explicitly 
only on the left hand side, evaluated at ri. This pro- 
vides many advantages in dealing with large field per- 
turbations, as will soon become apparent. In most cases 
we will consider perturbations about a uniform system 
with chemical potential p and density p{p) = p(r; p, [0]). 
When using this simplified notation p{p) should not be 
confused with p(r) = p(r; p B , [0]). Similarly, we will let 
p(p) denote the chemical potential of the uniform fluid 
as a function of density p. In a uniform system the direct 
correlation function c will take the simple form c(ri 2 ; p) , 
where r*i2 = |i"i — r 2 |. 

But how can we use Eq. to describe the density 
response to a large field perturbation such as the hard 
core field of interest here? This linear relation between 
a (possibly infinite) external field perturbation on the 
left hand side and the finite induced density change on 
the right must certainly fail for values of ri where the 
field is very large. Conversely, Eq. should be most 
accurate for those values of ri where the field is small — 
in particular where the field vanishes — and then through 
the integration over all r 2 it relates density changes in 
regions where the field vanishes to density changes in the 
regions where the field is nonzero. 

To treat large fields, we note that for any given ri we 
can locally impose the optimal condition that the field 
perturbation vanishes by introducing a shifted chemical 
potential 

A,* 1 = j* B - #n), (3) 

and a shifted external field 

^(r) = 0(r)-#n). (4) 

Since there is an arbitrary zero of energy and a constant 
external field acts like a shift of the chemical potential in 
the grand ensemble, we make no physical changes if we 
shift both functions by the same amount. In particular 

p(r; M B .M) = p(r; M* 1 ,^ 1 ])- _ 

The superscript ri in p Tl indicates a particular value 
of the chemical potential, which from Eq. J3Jl depends 
parametrically on ri through the local value of the field. 
When 4>(ri) vanishes, then p ri reduces to p B . We define 
p Tl , the hydrostatic density, by 

P^= P (t; f\[0]) = p(^). (5) 

Thus p ri is the density of the uniform fluid in zero field 
at the shifted chemical potential p ri ; equivalently p ri 
satisfies 

M(p ri ) = M ri =P B -Hn)- (6) 

Thus far, we have have merely introduced an equiva- 
lent (and apparently more complicated!) way of describ- 
ing the system in terms of a shifted field and a shifted 
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chemical potential. However this perspective immedi- 
ately suggests a very simple first approximation to the 
density response to a slowly varying external field. Since 
4> ri (r) by construction vanishes for r = ri , if (jf 1 (r) is suf- 
ficiently slowly varying, then the region around i"i within 
a correlation length is essentially in zero field. In that 
case the uniform hydrostatic density p Tl is clearly a good 
approximation to /?(i"i), the exact induced density at ri. 
Moreover, when the field is more rapidly varying, it is 
natural to introduce a second and even more accurate 
approximation to the density response. 

The hydrostatic density p Tl takes account only of the 
local value of the field at ri by a shift of the chemical po- 
tential. The HLR equation improves on this "local field" 
approximation by using linear response theory to deter- 
mine the density change from the hydrostatic density in- 
duced by nonlocal values of the shifted field <fi ri (r). Thus 
starting from the uniform density p Tl , we assume a linear 
response to the shifted field, replacing x~ ( r ij r 2! [p]) by 
X _1 (ri2;/9 ri ) in Eq. JIJ and setting S<fi(r) = <fi ri (r) and 
Sp(r 2 ) = p(r 2 ) - p Vl . Then the left side of Eq. JTJ van- 
ishes (giving the optimal linear response condition), and 
we have 



0= dv 2X -\r 12 :,p^)[p{v 2 )-p^ 



which can be rewritten exactly using Eq. J2J as 



P(ri 



Jdr 2 c(r 12 ; p ri )[p(r 2 ) 



(7) 



(8) 



This is our final result, which we refer to as the HLR 
equation. We view this as an integral equation relating 
the hydrostatic density p ri to the full density p(r) and 
solve it self-consistently for all r\. When 4>( T i) is known, 
we can immediately determine p ri at each ri from the 
local relation in Eq. ©, and then solve Eq. @ by it- 
eration for all ri to determine the full density response 
p(r). Conversely, for a given equilibrium density distri- 
bution p(r) we can use HLR equation to determine the 
associated field <f>{r) . This inverse solution of Eq. JBJ is 
particularly easy to carry out, since we can determine 
the local field at each ri separately, without iteration. 
Accurate results have been obtained for many test cases 
with strong repulsive or attractive fields 0, uj ■ 

This requires in particular expressions for p{p) and for 
the direct correlation function c{r\ 2 ;p) of the uniform 
hard sphere fluid. In this paper we will use the Percus- 
Yevick (PY) approximation for c(r 12 ; p) . This same 
function also arises from a self-consistent solution of the 
HLR equation, where the density response to a hard core 
field with R v = a (equivalent to fixing a solvent particle 
at the origin) is related to the uniform fluid pair corre- 
lation function. Thus this self-consistent use of the HLR 
equation provides a physically suggestive way of deriving 
the PY result for c{r 12 ; p) 0. The PY c(r 12 ; p) has a 
very simple analytical form and proves sufficiently accu- 
rate for our purposes here. Even better results can be 
found if one uses the very accurate expressions for the 



bulk c(r%2; p) and p(p) as given by the GMSA theory 
as inputs to the HLR equation. 



B. Relation to the PY approximation for a hard 
core solute 



A spherical cavity acts like a hard core external field 
that excludes the centers of all solvent molecules from 
the cavity region. We take the center of the cavity as 
the origin of our coordinate system, so that all distances 
are measured relative to the cavity center. Note that 
both the hydrostatic density p ri from Eq. (J5J and the full 
density response p(ri) from Eq. (|SJ) vanish whenever ri 
is located in the cavity. This exact "hard core condition" 
comes out naturally from the theory, and does not have 
to be imposed by hand as in the GFM or the GMSA 
approaches. 

To make contact with the PY approximation, re- 
call that the cavity-solvent direct correlation function 
C(i"i; p B , R v ) for this system exactly satisfies 



C(r i; p s ,i?„) = J drix- X {r l2 ;p B )\p{r2)-p B ]. 



(9) 



Thus C(ri) is the function that replaces — j38<t>(Yi) so 
that the linear response equation Q gives exact results 
when the full density change relative to the bulk is used 
on the right hand side. When this is compared to the 
HLR equation J7J for ri outside the cavity region (where 
p ri = p B and <j) = 0) we see that the HLR equation 
predicts that C(ri) vanishes. Thus for the HLR equation 
p{vi) vanishes inside the cavity region and C(ri) vanishes 
outside. This is the same as the PY approximation for 
the hard core cavity-solvent system 0, Ej ■ 

If R v is greater than a/2, with a the solvent hard core 
diameter, then an equivalent exclusion is achieved by re- 
placing the hard core external field by a hard core solute 
particle with (additive) diameter 



2R„ - a. 



(10) 



From this it follows that if the PY approximation for the 
bulk c or is used, the density p(r) predicted by the 
HLR equation is identical to that given by the PY equa- 
tion for the solute-solvent pair correlation function for a 
binary hard sphere mixture in the limit that the concen- 
tration of the solute species goes to zero 01 • Since an 
exact analytical solution of the PY equation for a binary 
HS mixture at arbitrary concentrations is known p"H , we 
can take advantage of these results when computing the 
excess grand free energy. 

This equality of solutions of the HLR equation and the 
PY mixture equation holds only for hard core cavity fields 
with radius R v > a/2 or a v > 0. As discussed below in 
Sec. for tiny cavities with R v < a/2 the HLR equa- 
tion can be solved directly and gives exact results for the 
density response if exact bulk correlation functions are 
used, and very accurate results when the PY approxima- 
tion for the bulk c is used. However, the corresponding 
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PY mixture solutions in this range of R v (arrived at for- 
mally by taking a v in Eq. (|10|) to be negative) are much 
less accurate. This inaccuracy arises from using the PY 
mixture solutions for negative a v . The direct solution of 
the PY cavity-solute equation for a tiny cavity, where a 
given approximation for the bulk c is used along with 
the PY approximation that C{y\) vanishes outside the 
cavity and p{y\) vanishes inside, gives the same accurate 
results as the HLR equation. However, for more general 
external fields, the HLR equation is quite distinct from 
the PY approximation, and is generally more accuratejit 
has given good results for a wide range of fields 0, la @ • 
This additional flexibility of the HLR equation will be 
required later in this paper when we discuss alternate 
density routes to the free energy. 



III. THERMODYNAMIC PATHWAYS TO THE 
FREE ENERGY 

In this section we discuss three different exact thermo- 
dynamic pathways for obtaining the excess free energy of 
inserting a cavity into a hard sphere fluid. The first two 
are well known, and the third describes a new density 
route that may have some computational advantages in 
other applications. We use the HLR equation to provide 
the needed structural information in all cases. We believe 
our calculation here represents the first use of a density 
route to obtain the excess free energy for this system. 
We then describe a simple but less accurate route to the 
free energy based on use of the partition function for the 
GFM. 



A. Compressibility route 

In this route the excess free energy is determined by 
varying the chemical potential of the system while the 
external field 0(r) producing the cavity with radius R v 
remains constant. In the grand canonical ensemble the 
average number of particles (N) is given by 



(11) 



where 0(/i, [4>]) is the grand free energy. We can then 
calculate the free energy difference between the final state 
of interest and the trivial ideal gas state of zero density 
with p = — oo and Q = by integration: 

V{H B , 14>}) = - dfi(N) = - / dn / drp(r;p, [</»]). 



(12) 

Then AQ V = £l{p B , [<f>]) - Q(p B , [0]), the desired excess 
grand free energy of the nonuniform fluid relative to the 
uniform bulk state, is given by 



As before, p(p) gives the density of the uniform hard 
sphere solvent fluid as a function of the chemical poten- 
tial. 

Since p{r;p, [<j)}) vanishes inside the cavity, Eq. fT3|l 
shows there is a term in the excess free energy propor- 
tional to the cavity volume v given by 

f^ B f pB dp 

v dpp(p)=v dp—p = vp B , (14) 

J-oo Jo dp 

on using the thermodynamic relation p(dp/ dp)r = 
(dp/dp)r- This exact leading order term for large v is 
determined when using the compressibility route so that 
p B is the uniform fluid pressure calculated by the com- 
pressibility route [l2| . 

The term in curly brackets in Eq. (|13|1 can be rewritten 
in a more convenient form for calculations by using the 
inverse relation to Eq. @ for a general chemical potential 
p: 

Pinw, [0]) - Pip) = J dr 2 x(ri 2 ] p)C(r 2 ; p{p), R v ), 

(15) 

where x( r i2! p) = P<>( r i — r 2) + P 2 [g{ r i2) — 1] is the usual 
linear response function of the uniform solvent fluid and 
g(r) is the radial distribution function. Substituting into 
Eq. i|13[) and carrying out the integration over ri , we have 
the formally exact result jl3l ] 

(3An v = -p f dp X {0,p) f dr 2 C{r 2 ;p(p),R v ) 



d P C(0;p, R v ). 



(16) 



AO, 



dp / dri{p(ri;/i, [<£]) - p( M )}. (13) 



Here x(0, p) is the k = value of the Fourier transform 
of x, with a similar definition for C(0; p, R v ). In the 
last equality we used the uniform fluid compressibility 
relation /?x(0, p) = dp(p)/dp to change variables to an 
integration over density. We will explicitly solve the HLR 
equation for R v < a/2 in Sec. Ivl below, and from the 
equivalence between the HLR equation and PY mixture 
equation for R v > a/2, we can use the exact solution of 
the PY mixture equation to obtain C at larger R v . Thus 
we can analytically carry out the integration in Eq. i|ltj|) 
for all R v . 



B. Virial route 

We now consider a different thermodynamic pathway, 
which was first used in scaled particle theory 0, llaj . 
Here we keep the chemical potential fixed at p B and vary 
the range of the external hard core field by a scaling 
parameter A, defining (f>\(r) = <fi(r/A). For the hard core 
cavity field of interest here, as A is varied from to 1 
the radius of the exclusion zone then varies from to R v . 
Since the density is generally related to the external field 
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in the grand ensemble by 
59. 



6(f>(r) 



(17) 



the desired free energy difference is given by integration: 



f3An v = dX drp x (r) 



dX ' 



(18) 



Here p\(r) = p(r;p B , {<j)\}). This formula is quite general 
and holds for any A-dependent potential that vanishes for 
A = 0. By exploiting special properties of the scaled hard 
core potential (the derivative of the Boltzmann factor of 
a hard core potential is a delta function) it is easy to 
show that Eq. I|18H can be exactly rewritten as 



/3AQ V = 4nRl l dXX 2 p\(XR. 



(19) 



Here p\{XR v ) is the contact density at the surface of the 
scaled exclusion zone with radius XR V . As in the com- 
pressibility route, we can analytically carry out the inte- 
gration in the virial route to obtain solvation free ener- 
gies for cavities for all R v . The equivalent PY solution 
for binary hard sphere mixtures is used for the contact 
densities for all Ai?„'s larger than a/2, while the explicit 
solution of the HLR equation is used for the XR V 's smaller 
than a/2. 



C. Density routes 



(or is strongly repulsive), then the region near A = 1 in 
the A-integration in Eq. lj2U|) must be treated carefully, 
since for r in the hard core region dp\(r)/dX is constant, 
while 4>\(j) must tend to infinity as A — > 1. Although 
the singularity in the potential is integrable (for a hard 
core potential the divergent term in (3(f>\ goes as — ln(l — 
A) and could be treated separately), in any case large 
contributions to the integral arise from a relatively small 
interval near A = 1. This could cause problems in a 
numerical integration. 

To illustrate the computational advantages and flexi- 
bility that different pathways can provide, we introduce 
here a new density route that removes this possible dif- 
ficulty. We consider a path that is linear in the square 
root of the density: 

p\" (r) ^ (p*)V* + \[pV\r) - ( P B )V% (22) 
where p^ 2 (r) = \J p\(r) , etc. For this pathway we have 



dpx{v) 
dX 



= 2^(r)&|a. 



(23) 



Both factors on the right side of Eq. (|23[) are easy to 
determine from Eq. (|22l) . The numerical integration in 
Eq. (|20|l can now be carried out straightforwardly since 

the p]/ 2 {r) factor in Eq. (|2*3|l will cause dp\(r)/dX to 
tend to zero exponentially fast wherever 4>\(r) becomes 
large. Results using this path are reported below. Other 
paths implementing this idea exist and we have not tried 
to make an optimal choice. 



In addition to these particular pathways, we can also 
imagine directly changing the equilibrium density from 
p B to p(r) over some convenient pathway specified by a 
coupling parameter A, while taking account of the asso- 
ciated changes in fl and </>(r). Integrating Eq. I|18|l by 
parts to make p\(r) explicitly the controlling variable, 
we have exactly 

0AQ v = J drp(T)<P(T)- J\x J dr<Mr)^l (20) 

Here <p\ (r) is the external field consistent with the spec- 
ified density profile, so that p(r; p B , [<f>\]) = P\(t). For 
a given density field p\(r), the HLR equation JSJ) can 
be solved inversely to obtain the associated hydrostatic 
density field p\. Using Eq. 10, <f)\ (r) at each r is locally 
related to p\ through p{p). Here we used the accurate 
Carnahan-Starling expression for p,(p). 

Most workers have considered a simple linear density 
path where 

px{v)=p B + X[p{v)-p B ]. (21) 

This has some theoretical advantages since dp\(r)/dX 
is independent of A and has been successfully used in 
numerical calculations of the surface tension of the liquid- 
vapor interface |£g. However, when cf>(r) has a hard core 



D. Gaussian Field Model 

Finally we consider an alternative approach, the Gaus- 
sian field model (GFM) 0, that for hard core fields has 
many common elements with the HLR method. The 
GFM describes density fluctuations in a uniform fluid 
with average density p by an effective quadratic Hamil- 
tonian 

H B = ^Jdr 1 J dv 2 6p(r 1 ) X - 1 (r 12 ;p)8p(r 2 ), (24) 

where 5p(r) = p(r)—p with p(r) the microscopic density. 
The partition function for a system in an external field 
<f)(v) = <f> (r) + 1 (r), with <fio a hard core field producing 
a cavity of radius R v and 4>\ a weaker perturbation, is 
then assumed to be given by [2j 



5, = Jm*) gw 

x exp[-(3H B +(3 J dr,5(r)0i(r)]. (25) 

The product of delta functions imposes the constraint 
that the density vanish inside the cavity. Inserting a 
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Fourier representation for the ^-functions and formally 
integrating p(r) from — oo to oo yields a Gaussian ap- 
proximation for the partition function, as discussed be- 
low. 

Moreover, using the same approximations, by function- 
ally differentiating "E v with respect to the field, one ob- 
tains the nonuniform singlet density in the GFM. In the 
case of a pure hard core field with <fii — 0, the density 
response to a cavity with radius R v is given by 



= p B - P 1 



' dr 2 dr 3 x{r 12 ;p B )x ln 1 (r 2l r 3 ) 

J V J V 



(26) 

The integrations are restricted to the cavity region, as in- 
dicated by the subscript v on the integral symbols. Here 
Xi n is the inverse of the restricted linear response func- 
tion Xin(ri2', H B ), which equals x( r i2', H B ) if both ri and 
r 2 are in the cavity region and equals zero otherwise. 
Thus Xin i s nonzero only inside the cavity and satisfies 



/ dr 2 x(ri2;p B )x t n( r 2,r 3 ) = <5(r a -r 3 ), 

J V 



(27) 



when both ri and r 3 are in the cavity region. Comparing 
Eq. I|26(l to Eq. i|I5|l . one can identify the cavity-solvent 
direct correlation function in the GFM as 



C{v 2 ) 



dr3X ln 1 (i-2,r 3 ) 



(28) 



By properties of Xin ■> t ne GFM C(r) vanishes outside 
the cavity region and p(r) in Eq. (|26|l vanishes inside. 
Thus the GFM gives exactly the same solution for the 
density response to a hard core external field as the PY 
or the HLR equations. (In the more general case where 
there is an additional perturbation potential (f>i , the var- 
ious approaches differ. The GFM can be shown to treat 
the softer tail using the mean spherical approximation, 
which is different from and generally less accurate than 
the hydrostatic shift used in the HLR equation.) 

Thus for cavities or hard core solutes all structurally 
based routes to the excess free energy will give the same 
results when using the GFM or the HLR equation. In ad- 
dition, the GFM partition function also provides a direct 
and very simple route to the free energy ;3j. However 
this route is inherently approximate because Eq. I|25(l is 
not really a free energy functional for the whole configura- 
tion space, but rather a restricted one describing only the 
space outside of the specified cavity region. This func- 
tional may legitimately describe subsequent small pertur- 
bations of 0i outside of the cavity, but it does not contain 
enough information about the functional dependence on 
the cavity volume in the first place. Moreover the approx- 
imations made in evaluating the GFM partition function 
do not build in the fact that in grand canonical ensemble, 
the thermodynamic properties should depend on [i — (f> 
rather than on p and </> individually. Thus it is also not 
consistent with the free energy prediction from the com- 
pressibility route, which integrates over states at differ- 
ent chemical potentials but with a fixed hard core always 
present. 



Evaluating the Gaussian integrals in Eq. (|25|l . the ex- 
cess grand free energy arising from a cavity with radius 
R v is given by 

PAfl v = -\og~ v /Z B 

= -lp B C(0;p B ,R v ) + log(det Xl n). (29) 

Here S„ denotes the partition function with no parti- 
cles in the cavity region and S B is the uniform bulk 
partition function. The term involving C in Eq. I|29|) 
arises from the Gaussian integration of p in Eq. i12">l) 
and Eq. (|28|l . When compared with the exact Eq. (|lfc>|> 
from the compressibility route, we see the GFM ef- 
fectively approximates C(0; p, R v ) for intermediate den- 
sity values by (p/p B ) C(0; p B , R v ). This free energy 
contribution has a form similar to a harmonic oscilla- 
tor with (p/p B ) C(0; p B , R v ) analogous to the restoring 
force. The second term is a result of the reduction of the 
configuration space. 

An alternate perspective considers the average proba- 
bility P V {N) of finding N particles in the volume v with 
radius R v . The probability of inserting a cavity is thus 
P„(0). Then a formally exact expression for the excess 
free energy /3AQ V is 



(3AQ V = ~ logP„(0) = - log N ^}°} — ■ 



(30) 



Here S^A] is the constrained partition function when 
N particles are in the specified volume. This formula 
has been successfully used in the information theory ap- 
proach developed by Hummer, Pratt and coworkers |l7j| . 

If the GFM is used to approximate the partition 
functions in Eq. (|30|l by replacing the product of 6- 
functions in Eq. I|25|l by the single average constraint 
S [J drp(r) — N], one arrives at a Gaussian approxima- 
tion [l8( for S„ [N] . This "discrete" approximation for 
(3AD, V based on this use of the GFM is 



(3An v 



l0g E«e-(«- s ) ! /^ 



where 



Xv = dri dr 2 x(|ri - r 2 |), 



(31) 



(32) 



and N = p B v. 

We determined the solvation free energy /3Afi„ for the 
GFM using both the continuum version, Eq. I|29|l . and 
the discrete version, Eq. I|31(l . For the uniform x w e used 
the PY result. To estimate logdet^m we expanded Xin 
in the volume v using two orthogonal basis functions. 
A single constant basis function was used in Ref. [Ls| . 
This is exact for tiny cavities with R v < a/2, as can be 
seen using Eq. 1|44J) below. We chose one basis function 
to be constant. The other was taken to be j$(R v r /it), 
the zeroth order spherical Bessel function with its first 
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node fixed at r = but made orthogonal to the first 
(constant) basis function. The second basis function is 
thus a linear combination of jo(R v r/Tr) and a constant. 
This was introduced to test the accuracy of the one ba- 
sis function approximation previously used and hopefully 
will give improved results for larger R v . 



IV. RESULTS FOR LARGER CAVITIES WITH 

Rv > a/2 

We now discuss the solvation free energies given by 
the various pathways for a cavity with R v > er/2, equiv- 
alent to a physically realizable hard core solute particle 
with diameter a v > 0. (Results for tiny cavities with 
R v < a/2 are discussed in Sec. Ivl below.) We use the 
simplest version of the theory, where the PY approxima- 
tion is used for the uniform fluid correlation functions. 
Fig. (|TJ give the solvation free energy /3AQ„ from the 
different pathways as a function of the packing fraction 
7] = np B a 3 /6 for R v /<r — 1, 1.5, 1.75, where the results 
can be compared to computer simulations °f Crooks 
and Chandler. Note that the volume of a spherical cav- 
ity with R v = 1.75cr is over 42 times greater than that of 
a solvent particle. For R v = a the results also give the 
excess chemical potential as a function of density for the 
uniform hard sphere fluid. 

As discussed above, we can obtain analytical expres- 
sions for AQ V for both the compressibility and virial 
routes. The compressibility route gives 



r?(-2 + lx] - II77 2 ) 
2(1 -V) 3 
18?j 2 (l + n) R v 2 



- log(l - v) 



18r] 3 R v 



(1 - ?y) 3 c 
and the virial route gives 



[iAVL v = 



r)(-2 + 7r? - hrj 2 ) 

I877 2 (1 - rj) R v 2 
(1 — 77) 3 a 2 



8 V (l + tj + V 2 ) R v 3 
(l-rj) 3 a 3 ' 

log(l - 77) 

877(1 + 77 - 2r] 2 ) R v 3 
(I-77) 3 a 3 



(33) 



(34) 



Results for the density route and for the GFM are com- 
puted numerically. 

In this range of R v there is good agreement except at 
the highest densities between the compressibility, virial 
and density routes, with best results overall arising from 
the compressibility route. The direct GFM predictions in 
Fig. (J2J from the partition function are less satisfactory. 
Both the discrete and the continuum versions of the GFM 
give results that approach zero incorrectly as p B — > 0, 
and the continuum values are consistently too large at 
high density while the discrete values are too small. The 
discrete version of the GFM uses a Gaussian approxi- 
mation for the constrained partition functions and gives 
less accurate results than could be obtained from a fit to 
accurate values of (N ) and (N 2 ) as in the information 
theory approach |l7j . 
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FIG. 1: The excess free energy predicted by the three ther- 
modynamic routes for cavity radii R v = a, 1.5a and 1.75a, 
compared with simulation data. 77 is the packing fraction and 
is equal to irp B a 3 /6. 
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FIG. 2: The excess free energy predictions by both the dis- 
crete and the continuum versions of the GFM, plotted for the 
same R v values and on the same scale as in Fig. 



As R v — * 00, the surface of the cavity approaches that 
of a planar wall. As shown in Eq. (|14f) there is a diverging 
term in the excess free energy given by the cavity volume 
v = 47ri? 3 /3 times the bulk pressure p B , and the more 
interesting quantity to calculate is the surface term j v , 
given by 



Plv = 



f3An v ~ fip B v 
AttR 2 ■ 



(35) 



The surface tension of the planar wall is then 700. In 
the present case, both the compressibility and the virial 
routes give analytical expressions for /3AQ, V which depend 
oni?„ as a polynomial: ao+aiR v /a+a2Rl/<J 2 +a3R 3 ,/a 3 . 
The coefficient 03 thus gives another route to the bulk 
pressure on taking the wall limit. The p B used in Eq. I|35|) 
has to agree with the prediction from the 0,3 so that 
is finite. As discussed earlier, the compressibility (3Afl v 
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FIG. 3: Surface tension predicted by the thermodynamic 
routes compared with the simulation fitting formula. 



yields the same bulk pressure as given by the accurate 
uniform fluid PY compressibility equation of state. The 
virial route does not automatically build in this consis- 
tency, and the pressure predicted from the coefficient 03 
is less accurate than the uniform fluid PY virial equation 
of state. 

The coefficient of the quadratic term then gives the 
surface tension /3 7oo = 0-2 /47tct 2 . Using the compressibil- 
ity route we find 



- 47r/3 7oo cr^ = 
while the virial route gives 



47T/3 7o 



2 _ 18r7 2 (l + r?) 



(1 - nf 



18?? 2 (1 -rj) 
(l-ry) 3 ' 



(36) 



(37) 



The 7oo obtained by the compressibility route coincides 
with that given by scaled particle theory 0, 0, |2(j ■ We 
obtained 700 numerically for the density route. 

These results can be compared to the quasi-exact for- 
mula pn 



18?7 2 (1 



44 

35*? 



h 2 ) 



(38) 



which fits simulation data [20| and imposes the known 
first and second surface virial coefficients j^. As shown 
in Fig. ((2J) the compressibility and density routes give 
excellent results, while the virial route is much less sat- 
isfactory. 

This can be understood since the virial route uses only 
the contact densities at the fixed bulk density. The HLR 
equation is least accurate for the contact density at high 
bulk density and for large R v , while the density response 
away from the solute is more accurate. On the other 
hand, the compressibility and density routes make use 
of the density response at all distances and over a range 
of densities from low density to the final p(r) where the 
HLR equation is more accurate. 



V. A SPECIAL REGIME: TINY CAVITIES 
A. Exact results 

The density response to a tiny cavity with R v < a/2 
is especially simple, since the center of at most one hard 
core solvent particle can lie anywhere within such a region 
[Til Hil . This fact allows one to determine exactly both 
the density response to a tiny cavity and the solvation 
free energy. Here we will compare these exact results to 
the predictions of the HLR and GFM approaches 

As in Eq. 13U|) . the solvation free energy is directly 
related to the average probability that no particle are in 
the cavity region 



Pv(0), 



(39) 



where P V (N) is the probability of finding TV particles 
simultaneously in the region with volume v = AttR^/3. 
When R v < a/2 , the region can hold no more than one 
solvent particle, so that 



and 



P V (0) + P V (1) = 1, 



p B v = (N) V =P V {1). 



(40) 



(41) 



Substituting these into Eq. I|39|l we thus find the exact 
result |3 Eg 



(3An v = -\og(l-p B v). 



(42) 



This argument can be extended to show that the exact 
density response to a tiny cavity is |23|: 



P(r) 



p B {l-f v dr>p B g(\r> 
(l~p B v) 



(43) 



for r outside v, with p(r) = for r inside. Here g{r) 
is the exact radial distribution function for the uniform 
solvent fluid. Note that the contact density p(R v ) — 
p B /(l — p B v) is exactly determined independent of the 
details of g(r), since the corresponding g(\r' — r|) in Eq. 
(|43H vanishes for all r' inside v. This result is valid as 
long as the inserted region v can hold no more than one 
solvent particle, so Eq. (|43|l also holds for solvents with 
a hard core pair potential plus a softer tail. 



B. Structural predictions of the GFM and HLR 
methods 

Now let us examine the GFM result in Eq. (|26[1 for 
the case of a density response to a tiny cavity. Since 
g(\ri — i"2 1 ) = when both 17 and r 2 are in v, then 
X™(ri,r 2 ) = p B S(r 1 - r 2 ) - (p B f inside v. It is easy to 
see from Eq. <(T7I) that the inverse function xin t nen has 
the simple form |24| 



(ri,r 2 ) 



5(i7 - r 2 ) 



1 



p B v 



(44) 
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FIG. 4: The compressibility and virial routes are exact for 
the tiny cavity regime. The density route is plotted along 
with the GFM results to compare with the exact free energy 
predictions for R v = 0.3<r and 0.5a. 

When Eq. (|4*4"|) is inserted into Eq. j^H} to obtain the 
GFM density response to a tiny cavity, we recover the 
exact expression for p(r) given in Eq. I|43|l . provided that 
the exact uniform fluid g(r) or c(r) is used. If approxi- 
mate (say PY) results are used to describe the uniform 
fluid response functions then strictly speaking the GFM 
and HLR predictions for p(r) will not be exact for all r. 
However, the contact density p(R v ) is exact in any case, 
since, as noted earlier, this requires only that the ap- 
proximate g(r) vanish inside the cavity region. Because 
of the equivalence between the structural predictions of 
the GFM and the HLR equation, these same conclusions 
hold for the HLR equation. In particular, the density re- 
sponse outside a tiny cavity is exactly described by linear 
response theory about the uniform bulk system. 

C. Solvation free energy from thermodynamic 
pathways 

Moreover, these theories give the exact result of Eq. 
(l4"2|l for the solvation free energy (3AQ V , independent of 
possible errors in g(r), for all structurally based thermo- 
dynamic pathways that use only tiny hard core fields. In 
particular, the virial route in Eq. i|19|l gives exact results 
for /3Af2„ since this requires only the (exact) contact val- 
ues p\(XR v ). The compressibility route as written in Eq. 
(|16|l requires C(ri ) , which from Eq. (|28() depends only on 
the exact \in m -^q- (|44fl . Both results require only that 
the bulk g(r) vanish for r < a, and are unaffected by any 
errors at larger r. This is in accord with our general sup- 
position that particular thermodynamic pathways can be 
relatively insensitive to errors in the structural theory. 

However the choice of pathway is important. Thus the 
density routes do not give pAQ v exactly even in the tiny 
cavity regime. This is because as p\(r) is varied, the 
corresponding <p\ (r) in general is not a pure hard core 



field and spreads outside the cavity region. Neither the 
GFM nor the HLR theories can treat these softer and 
longer-ranged fields exactly even if exact uniform fluid 
correlation functions are used. 



D. Solvation free energy from the GFM partition 
function 

The j3AVt v obtained directly by taking the logarithm 
of the partition function in Eq. I|29[l can be expressed 
analytically for tiny cavities as 

(3AQ V = \ log (1 - p B v) + -tr^r • (45) 
2 2 1 — p a v 

As v — » 0, this goes as (p s u) 2 /4, while the exact result 
from Eq. (|42|l goes as p B v. This quadratic term arises 
from the assumption of Gaussian fluctuations, which 
breaks down in this limit. The alternate discrete version 
from Eq. i|3T|) reduces to 

e _iV/[2(l-iV)] 

(3AU V = - log ^_-_— — 7pT , (46) 

where N = p B v. This has the peculiar behavior as v — > 
that all derivatives vanish, and so is significantly in error 
in this regime. See Fig. 21 for comparison with the exact 
answer. 



VI. CONCLUSION 

We have discussed several different thermodynamic 
routes that can be used to determine the solvation free 
energy for inserting both small and large cavities into a 
hard sphere fluid. Generally accurate results are found 
by using the HLR equation to relate the densities and as- 
sociated fields over the intermediate states of the differ- 
ent pathways. We also considered the GFM and showed 
that it gives results equivalent to the HLR equation for 
the density response induced by a rigid cavity. However 
the GFM cannot describe the softer external potentials 
and intermediate densities needed for the density routes 
and for more general thermodynamic pathways. Direct 
use of the approximate partition function of the GFM to 
determine the solvation the free energy of a cavity also 
gives less accurate results. 

Best results using the HLR equation for the solvation 
free energy of a cavity are found from the compressibility 
and density routes. This can be understood since most 
states along these routes require the density response at 
intermediate densities and distances away from the cav- 
ity where the HLR equation is most accurate. The HLR 
equation can also be used for more general solutes with 
different shapes or longer ranged attractive interactions 
and in applications where other pathways may be more 
useful. Combined with an appropriate pathway it repre- 
sents a versatile and computationally efficient method for 
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determining both the structure and the thermodynamics 
of nonuniform fluids. 
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Here p B G(R) is the contact value of the density response 
to a cavity of size R. In the case considered in this pa- 
per, the 7oo from both the compressibility and the virial 
routes can be determined analytically. In general, when 
no analytical expressions are available, one may need to 
carry out these integrations numerically. 



APPENDIX 



In the "wall" limit where R v — > oo, the surface tension 
given by the compressibility route is determined from 



7oo 



dp ^T dz[p w (z)-p], (A.l) 
o °P Jo 



with p w (z) = p(z + Rv). The virial route gives 

/•OO 

loo = k B Tp B dR[G{R) - G(oo)], (A.2) 
Jo 



ft 



There also exists the exact virial sum rule p B G (oo) = 
B for a planar wall immersed in the hard sphere fluid 
|25j . However, this focuses on the structure at contact for 
the planar wall, which is the region where the HLR equa- 
tion is least accurate. Thus this thermodynamic pathway 
for the bulk pressure gives relatively poor results. 
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